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Abstract: We solve the CCFM equation numerically in the presence of a boundary condi- 
tion which effectively incorporates the non-linear dynamics. We retain the full dependence 
of the unintegrated gluon distribution on the coherence scale, and extract the saturation 
momentum. The resulting saturation scale is a function of both rapidity and the coherence 
momentum. In Deep Inelastic Scattering this will lead to a dependence of the satura- 



tion scale on the photon virtuality in addition to the usual xbj dependence. At asymptotic 
energies the interplay between the perturbative non-linear physics, and that of the QCD co- 
herence, leads to an interesting and novel dynamics where the saturation momentum itself 
eventually saturates. We also investigate various implementations of the "non-Sudakov" 
form factor. It is shown that the non-linear dynamics leads to almost identical results for 
different form factors. Finally, different choices of the scale of the running coupling are 
analyzed and implications for the phenomenology are discussed. 
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1. Introduction 

The linear QCD evolution equations at small-x, such as BFKL [1] or CCFM [2-6], lead 
to a very rapid growth of the gluon density. Once the density is too large however, one 
needs to take into account the non-linear effects. The BFKL evolution is modified in 
the high density limit by the so-called Color Glass Condensate (CGC) [7, 8] (for a review 
see [9]) formalism where the strong color fields are treated semi-classically, and where 
their presence leads to the saturation of the gluon emission rate 1 . The BFKL evolution 
equation is then generalized to either the Balitsky-Kovchegov (BK) [10,11] equation or to 
the more complete JIMWLK [12, 13] equation. The energy dependent momentum scale 
characterizing the transition from the weak to strong field regimes is called the saturation 
momentum and is denoted Q s (x). 

What makes CCFM different from the standard BFKL-type picture of the small- £ 
regime is the property of the so-called coherence of emissions. This is quantum mechanical 



1 This implies that the gluon density keeps growing also in the saturated regime but the growth is 
logarithmic in x as opposed to power-like. 



in origin and comes from the interference in the emission of gluons from different partons. 
This leads to what is termed angular ordering which, as a fundamental feature of QCD, 
is present in both time-like and space-like evolution as well as at both small and large 
x. In CCFM the effects of angular ordering are encoded in an additional momentum 
scale, usually denoted p, so that the unintegrated gluon density A now depends on three 
variables: Bjorken x, the momentum k, and the additional momentum p. We will denote 
the unintegrated gluon density by A(x, k,p). In comparison, the BFKL gluon density, and 
similarly the density in the CGC formalism, depends only on x and k, that is A = A{x, k). 

The saturation scale Q s is the momentum at which the gluons in the hadron start to 
overlap with each other. This leads to the increased rate for the recombination of gluons. 
Q s can be defined as the momentum for which A is a given constant, say of the order 
of 1. This is a natural definition since the gluon density can be defined as a phase space 
occupation number 2 of gluons per unit rapidity with given transverse momentum k. This 
automatically implies that Q s is x dependent. However, it is then clear that in CCFM 
the same definition immediately leads to an extra scale dependence of Q s , since now the 
requirement of constant occupation number, A(x,Q s ,p) = const, leads to Q s = Q s (x,p). 
One of our main objectives in this paper is to study this novel scale dependence of Q s by 
numerically solving the CCFM equation in the presence of a boundary condition enforcing 
non-linear corrections. 

The additional momentum variable p in CCFM is determined by the kinematics of 
the hard scattering. More precisely, the hard scattering determines the maximal angle 
of the allowed emissions in accordance with angular ordering. Conventionally one uses 
the squared angle £ = q±/qf as the angular variable, and for example in a DIS process 
with n gluon emissions, angular ordering implies that £1 < £2 < • • • < £n < £ where £ 
is the maximal angle determined by the hard scattering. The transverse momentum p is 
then defined via £ = p 2 /(x^s) where x n is the energy fraction of the gluon entering the 
hard scattering, and s is the cms energy. The variable p is related to Q in DIS, roughly 
£ — Q 2 /{x 2 s). This means that the saturation scale from CCFM in DIS kinematics will 
depend on the photon virtuality Q in addition to the Bjorken x. In the usual picture 
of the small- x physics, however, Q s is an intrinsic property of the target hadron, and is 
independent of the probe and its virtuality. The difference now is that the gluon occupation 
number in the hadron depends additionally on the probe. That Q s is not a property of the 
target alone but also depends on the probe and its virtuality is certainly what one would 
generally expect from quantum mechanics. In the end Q s is determined by the scattering 
process, to which both the probe and the target enters, and therefore we should expect Q s 
to depend on both. In that sense the extra scale dependence of Q s is most natural. 

Unlike BFKL, the proper non-linear generalization of CCFM is not know yet, and we 
shall not present it here either. In [14,15], however, a method was proposed which effectively 
implements saturation and unitarity corrections within a generic linear evolution equation, 
and the method was in particular applied to BFKL and CCFM. The idea is simply based 
on enforcing a boundary condition in the linear evolution which prevents the gluon phase- 

2 The exact relation between A and the occupation number depends on the precise definition of the 
former. 



space occupation numbers to grow too large, and the procedure is an extension of a strategy 
originally introduced in connection with analytic studies of BFKL evolution in the presence 
of saturation effects [16, 17]. The method was demonstrated to work very well for BFKL 
in [14] where the solutions of BFKL in the presence of the saturation boundary condition 
were compared to the solutions of BK, and it was shown that the method reproduces 
successfully the results from BK for any x and for both fixed and running a s . 

The method was then applied to CCFM in [15] and Q s was extracted. However, in 
order to simplify the numerical procedure, only the case k <C p was studied in which case 
one can neglect the dependence on p so that Q s depends only on x. This is the region that 
leads to another formulation of CCFM called the Linked Dipole Chain (LDC) model [18] 
which has some advantages compared to CCFM such as projectile-target symmetry. In 
this paper, we will keep the full p dependence of the evolution, and we will show that this 
dependence, which as mentioned above is related to coherence in the emissions, combined 
with the physics of saturation gives rise to a very interesting and novel dynamics which 
eventually leads to the saturation of the saturation momentum Q s itself. This happens 
when the solution to A "stalls" in some region of k, meaning that the evolution in Y is 
almost completely stopped. This behavior is in contrast to the standard picture of the 
non- linear dynamics where Q s grows indefinitely with Y. 

Generally we shall see that the shape of Q s and that of A close to the saturation 
regime is in CCFM rather different than in BFKL. Also some characteristic features of 
the standard small- a; evolution such as geometric scaling [19] is modified in CCFM. On 
the practical level, however, it seems that the new behavior of Q s , whereby it eventually 
saturates as a function of Y for given p, sets in very late, especially with a running coupling, 
so that it seems it would be be very hard to observe at collider experiments. On the other 
hand it will still be important to take into account the variation of Q s with p also for 
phenomenological values of Y, especially when p is small. 

In this study, we restrict ourselves to the version of CCFM containing only the "hard" 
emissions associated with the 1/z pole in the gluon splitting function. In addition to 
these, in CCFM there are also the "soft" emissions which are associated with the opposite 
1/(1 — z) pole 3 . Within the accuracy of the formulation of CCFM however, these emissions 
are exactly "probability conserving", meaning that they are exactly compensated by the 
associated virtual corrections encoded in the so-called "Sudakov" form factor. Thus at 
least formally, it should make no difference for the solution of A whether one includes 
the soft emissions or not. The difference would be only in the exclusive final state, a 
correct description of which requires the soft emissions. In practice, however, such a formal 
statement carries little significance, and the inclusion of soft emissions can have substantial 
effects on A. Thus it is our aim to also include the soft emissions in our procedure. This 
is, however, somewhat challenging numerically due to the fact that the real and virtual 
terms are separately divergent, and so one must introduce a cut-off which complicates the 



3 In addition to the singular pieces of the splitting function, the non-singular pieces are usually included 
in phenomenological applications, such as in the CASCADE Monte Carlo [20]. While these terms may 
obviously be important in practice it should be noted that their inclusion is somewhat ad hoc since they 
were not properly included in the original formulation of CCFM. 



numerical procedure. Therefore the inclusion of the potentially important soft emissions 
and the associated Sudakov form factor is postponed to a future work. 

The \jz emissions are accompanied by virtual corrections which are encoded in the 
so-called " non-Sudakov" form factor. One can in the literature find different expressions 
for this form factor. In this paper we shall study both of the common expressions for the 
form factor. We shall see that, generally, the different choices of the non-Sudakov form 
factors lead to rather different results in the linear evolution, which is due to the reason 
that the differences which occur in the softer region are enhanced unless there is some 
mechanism to cut off the growth there. As we will demonstrate, once non-linear evolution 
is included via the saturation boundary condition, then the differences are to a large degree 
removed. The small differences remaining when Y is extended to very large values can can 
be removed by scale factors independent of Y and p. Thus saturation introduces a certain 
universality in the evolution as was already observed in [15]. 

Related to this, it is important to check the sensitivity of the results to different 
prescriptions of the boundary condition. Generally, one can consider different boundary 
conditions or, for a given type of boundary, different values of the boundary parameters. 
In [15] it was reported that the different choices lead to different normalizations of the 
solutions, without altering the shape or the Y dependence. Since we here keep the p 
dependence, the question is more subtle. We will show that the differences are small and 
that they can generally be removed by Y and p independent scale factors. Some differences 
can remain in the region where Q s starts to saturate (the region where A "stalls"), but 
what we find is that this behavior itself, i.e the saturation of Q s , is independent of the 
precise application of the boundary. What changes is simply at what Y, for a given p, 
this behavior sets in. To know the exact details of the behavior in the saturated regime 
one would need to fully formulate the proper non-linear generalization of CCFM which is 
beyond the scope of this work. It seems, however, rather reasonable to expect that this 
interesting dynamics where, as a result of the combination of saturation and coherence, 
Q s eventually saturates will remain true also once the proper non-linear generalization of 
CCFM is known. 

The structure of the paper is the following: in the next section we present a brief 
overview of the CCFM formalism and show the numerical results for the solution in the 
case when the full dependence on the scale p is included. We also discuss different versions 
of the running coupling prescriptions. In section 3 we provide semi-analytical results for 
the case of CCFM with the boundary condition. In particular, we discuss the dependence 
of the saturation scale on the rapidity and p. In section 4 we briefly describe the method of 
implementation of the saturation boundary. In section 5 we present the complete numer- 
ical analysis of CCFM in the presence of the saturation boundary. We analyze different 
prescriptions for the non-Sudakov form factors and the running coupling, and we extract 
the rapidity and p dependence of the saturation scale. Finally, in section 6 we state our 
conclusions and discuss the relevance of the results for the phenomenology. 
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Figure 1: The kinematics of the gluon ladder. The energy fraction of the real s-channel and 
virtual i-channel gluons are denoted by j/j and x% respectively while their respective transverse 
momenta are denoted by qi and fcj. 

2. Overview of CCFM 



Our aim in this section is to provide a brief overview of CCFM, focusing on the points which 
shall be relevant for our analysis. It is possible to find different versions of the non-Sudakov 
form factor in the literature but the origin of this form factor is hardly ever discussed. It 
will therefore be important to recall the derivation of this form factor and understand the 
motivation behind the formulas. Moreover, we shall in addition to the standard choice of 
the scale of the running coupling also consider a different choice. The most comprehensive 
overview is to be found in the original papers in [2-6], and a simplified but rather detailed 
discussion on CCFM is also offered in [15]. 

The gluon ladder generated by the CCFM evolution is represented in figure ||. In 
CCFM the successive emissions along the ladder are ordered in their angle. As mentioned 
in the introduction, it is customary to define £ = Q±/qf as the angular variable (with the 
direction of the incoming proton defining the longitudinal direction) , and angular ordering 
implies that £j > £j-i. Additionally there is a maximum angle £, set by the kinematics of 
the hard scattering, which limits all emissions. 

In CCFM, the emissions are divided into hard and soft emissions, and the importance 
of this classification is that it offers significant simplifications in the formulation. The 
precise definition is that the hard emissions are those emissions ordered in both energy 
and angle, while for each soft emission there exists another emission with bigger angle and 
larger energy fraction. It follows that, the hard emissions are those associated with the 



1/z pole in the gluon splitting function while the soft emissions are associated with the 
opposite 1/(1 — z) pole. Following the above definition one can associate with each hard 
emission k a cluster of soft emissions Sk which satisfy £k-i < £ < £fc an d y < Uk where £k-i 
is the angle of the previous hard emission. The label "soft" has a two- fold meaning. Firstly 
the fact that each of these emissions takes away very little energy, because of the 1/(1 — z) 
pole, and secondly that all soft emissions in the set Sk satisfy q < %. The classification of 
the emissions as being hard or soft is also related to another classification of the ladder in 
terms of so-called k±— changing and k± -conserving emissions respectively. What is meant 
by this is that there are certain real emissions which leave the virtual propagator momenta 
unchanged in an approximate sense. This latter classification is crucial for deriving the 
'standard' version of CCFM which is used in all practical applications. 

2.1 Integral equations and the virtual form factors 

The CCFM integral equation, including both the soft and the hard emissions, is given by 

A(x,k,p) = a s / dz ^2 9 (p-zq) A s (p,zq)(—^ — — - \- — J A (-,k',q 

(2.1) 

The momentum variable p is defined via £ = p 2 /( x n s )> and A/ = |fc + (l — z)q\. The momen- 
tum q is the rescaled momentum of the real gluon, and is related to the true momentum q 
by q = q/{l — z). In line with all the approximations made in order to derive this equation 
one can for the 1/z pole set q = q (which is what we have already done in A ns (z, k, q)). 
Here A s is the Sudakov form factor which screens the singularity of the 1/(1 — z) pole, 
while A ns is the " non-Sudakov" form factor which depends also on the momenta k and 
it corresponds in BFKL to the gluon Regge factor. The rescaled coupling constant a s is 
defined as a s N c /ir. 

In figure |2| we represent pictorially the components of equation (2.1). Here the horizon- 



tal axis represents the inverse energy fraction of the gluons, the further to the right a gluon 
is located the less energy it has, while the vertical axis represents the transverse momenta, 
the higher up a emission is located the larger its transverse momentum. The maximal 
angle ^ is represented by drawing a line of constant angle from the point (x n ,p); this is the 
diagonal line going through the black dot representing p (p itself is not the momentum of 
any gluon (real or virtual) but we still denote its location by a dot in the figure). In this 
particular case, k < p. The "Sudakov" form factor A s (p,zq) is then given by exp(— a s 5) 
where S is the phase space area which is the diagonally shaded region in the figure. This is 
the region where further emissions would have been possible. Similarly the "non-Sudakov" 
form factor A ns is given by exp(— a s A) where A is the area of the vertically shaded region. 
In the figure we also show the next backward step in the evolution which is obtained from 
(|2.l| ) when A in the RHS is expanded iteratively. As we can see from the figure, the region 
S contributing to A s diverges when extended down indefinitely. Thus this divergence is cut 
by a horizontal line in the figure, which means a soft cut-off in the transverse momentum. 
Notice that this divergence is the one arising from the 1/(1 — z) pole because as z — > 1, 
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Figure 2: Graphical representation of equation (2.1). The horizontal axis denotes the inverse 
energy fraction while the vertical axis denotes the transverse momentum, all in logarithmic scale. 
The shaded regions show the regions over which the form factors receive contributions. The real 
emission q is represented by a dot, while the virtual gluons k and kl are represented by crosses. 
The position of p is also indicated. The diagonal lines are lines of constant £. 



the gluon carries infinitely small energy fraction, and in the figure it would then be located 
precisely in the region where the diagonally shaded region is extended down indefinitely 
(the fact that this region is bounded by the diagonal lines is due to angular ordering). 

The soft emissions are such that they are exactly compensated by the factors A s . 
This means that the regions over which the Sudakov form factors receive contribution 
are precisely the regions to which the soft emissions are confined, i.e. the clusters Sk 
mentioned above. Summation over all possible soft emissions in each Sk leads to a factor 
exp(+a s Sfc). Each such factor then multiplies to unity with the corresponding Sudakov 
form factor exp(— a s Sk). One is consequently left with the hard emissions only, and the 
integral equation reads 



A(x,k,p) 



a\ 



dz f d 2 q 

irq 2 



d{p - zq) A ns (k, z,q)A[ -, k', q 



(2-2) 



This equation is thus formally equivalent to (|2.1|). In practice, however, one can expect 
to find different solutions. The numerical implementation of the soft emissions faces the 
difficulty that one must cut the z — >• 1 singularity by a momentum cut-off as in figure |2|. 
This introduces a numerical uncertainty which complicates the procedure. Therefore, even 
though their inclusion is important and interesting, we will postpone their treatment to a 
future work and in the present concentrate only on equation (|2.2|). 

As mentioned above, the non-Sudakov form factor A ns can be written as exp(— a s A), 
where in figure ||, A is the vertically dashed region bounded by k and the angle of q, and 
the energies of q and k. In CCFM, the form factors A ns and A s are derived out of the so- 
called eikonal, S e n~, and non-eikonal, S ne , form factors which sum up to all orders the virtual 
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Figure 3: Illustration of the regions contributing to the definition of the non-Sudakov form factor 
A ns . In this case k > q. 



corrections associated with the eikonal and non-eikonal emission vertices respectively [2-6] . 
As shown in [15], for each emission k the product of the associated form factors S e ik{k) and 
S ne (k) equals the BFKL form factor A BFKL (k), that is S eik {k) ■ S ne {k) = A BFKL {k). The 
definition of A ns for a given emission q is that it is the product between S ne and a portion 
of the eikonal form factor which covers a region bounded from below by the angle of q and 
from above by the maximal angle £ (this is the region left over after multiplication with a 
soft factor coming from the resummed real emissions). In figure ^ this is the sum of the 
shaded regions B\ and C\ . The non-eikonal form factor S ne itself is such that it covers a 
region bounded from below by k while from above it is again bounded by £, this is region 
B\ in the same figure. What is important however is that the sign in the exponential factor 
is reversed, so that it is positive. 4 The non-Sudakov is then given by 

A ns = exp(-a s (B 1 + C{)) ■ exp(+a s Bi) = exp(-a s Ci) . (2.3) 

In this case, however, we assumed that k > q, see figure pJL Now, if k < q we can have 
two different situations, depending on whether or not k is entirely below the diagonal line 
through q in the figure, that is whether k < zq or k > zq. In these both cases we define 
a new region A\ shown in figure Q. This is a region bounded from below by k and from 
above by the angle of q. When k > q as in figure || of course A\ = 0. 

On the left plot of figure || we show the case when k > zq, and the definition of A ns 
then yields 



A ns = exp(-a s (B 2 + C 2 )) ■ exp(+a s (^ 2 + B 2 )) = exp(-a s (C 2 - A 2 )) . 



(2.4) 



4 This actually comes from the fact that S ne contains virtual corrections which are generated by the 
product between the non-eikonal current J ne and the eikonal current J e ik, while S e ik is generated by the 
square of the eikonal current J^ ik . This leads to a sign difference between the form factors, for details 

see [2-6]. 
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Figure 4: Same representation as in figure H but in this case with zq < k < q (left) , and k < zq 
(right). 



If A2 is bigger than C2 then we see that A ns is larger than unity. In case k < zq we instead 
have (in this case C 3 = 0) 



A ns = exp(-a s B 3 ) ■ exp(+a s (^ 3 + B 3 )) = exp(+a s A 3 ) > 1 . 



(2.5) 



Thus A ns exceeds unity when A{ > d, and this actually happens when k 2 < zq 2 . All the 
results above give 



A ns = exp -a t 



1 dz , r k 
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(2.6) 



Therefore if the kinematical constraint k 2 > zq 2 is assumed to be valid, A ns is guaranteed 
to be bounded by unity 5 . A different version of A ns was presented in [21] in which case 
the regions Ai are simply neglected. Then one always has 



A ns = exp(-a s Ci), 



where 



1 k 2 
C\ = In - In — ^ 

z zq z 



Co = In 2 



C* = 



(2.7) 



(2.8) 



zq 



Since Ai is neglected, we are thus guaranteed that A ns < 1. Thus in a sense this has a 
similar effect to assuming the kinematical constraint to hold, but one should still keep in 
mind that the effect is not the same because one has no apparent reason to simply ne- 
glect the regions Ai (this amounts to neglecting soft emissions above k, see also discussion 
below). What the kinematical constraint ensures is that Cj > Ai. In this paper, we will 



5 In that case A ns , despite its name, works as a Sudakov and it can be used to cancel the fc^ -conserving 
hard emissions. One is then left with a much simplified formula. This was first understood and used in [18], 
and it was used in the analysis in [15] as well. 



not include the kinematical constraint (it was included in [15]), which must be included in 
the real emissions as well, but we shall rather study the differences between using the two 
different non-Sudakov form factors, ( |2.6D and (|27~ 



2.2 Scale of the running coupling 

Above we assumed a fixed coupling which is consistent with formulation of CCFM which 
was only performed at LO where the coupling is fixed. Although the effects of the running 
coupling are formally of next to leading order, it is nevertheless for any phenomenological 
application extremely important to take into account the running of the coupling. 

The standard procedure is to let the coupling run with k in the 1/z pole and in A ns 
which is similar to what is usually done in LO BFKL. When both hard and soft emissions 
are included then q is usually chosen as the relevant scale for a s in the 1/(1 — z) pole and 
in A s (this is for example done in CASCADE as well as in the SMALLX program [22]). 
Notice that A s multiplies also the 1/z pole so with this choice one is for the hard emissions 
using q as scale for the Sudakov form factor while k is used for the non-Sudakov and the 
real emissions. 

As k is the larger scale in A ns it may seem natural that it is chosen as the scale for 
the running coupling. However, once a scale is chosen for a certain class of emissions, 
say the hard emissions, then that scale should be consistently applied everywhere because 
otherwise some cancelations which are important for the consistency of the formalism may 
be violated. Also, even though k seems to be the larger scale in A ns we have seen above 
that actually in the derivation of A ns , k first enters as the lower scale in S ne . 

In [22] it is argued that the choice to have q and k as scales simultaneously can be 
derived from having the single choice (1 — z)-y/|fc 2 | (where k 2 is now the four momentum 
squared), since this choice reduces to k and q when z — > and z — >• 1 respectively. One 
should, however, again remember that both A ns and A s are derived out of S e ik and S ne , 
and while the four momentum of k could be argued as scale in the latter, the former resums 
the graphs where virtual gluons are emitted and absorbed in between the real gluons so it 
is hard to see why the scale of k should also be used here. We find it more plausible that 
instead the momenta q of the virtual gluons are used in S e n-- But if this choice is made, 
then to derive A ns and A s in the first place, one has to choose q as scale in S ne too, since 
otherwise necessary cancelations between S e ik and S ne would not work. Consequently this 
leads to using q in both A ns and A s , and therefore also in all the real emissions. 

If one uses k as scale in A ns then obviously the formula in ( |2.6[ ), and similarly the 
corresponding one in ( |2.7D is the same as in the fixed coupling case. If on the other hand 
one uses the virtual momentum q' as scale for a s , then the formula is changed. We use q' 
as scale in ( |2.7| ) and derive a new formula which is however somewhat lengthy so we do 
not present it explicitly. We shall use the one loop result for the running coupling, frozen 
below some scale ko, 

a s {q 2 ; k 2 ) = - 2 -, q = max(g, k ) ■ (2.9) 
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Thus for the running coupling case we shall consider the following alternatives: One in 



which we use k as scale and A ns given in (2.6), one in which we use k as scale and A 



in (|2.7| ), and the last case where q is chosen as scale in a s and the corresponding modified 
A ns is used. 

We shall note here that, in the BFKL equation the coupling also starts to run at the 
next-to- leading logarithmic level (in lnl/x) only [23-28]. The different choices of scales 
formally lead to the differences at next-to-next-to leading level. They are however not 
negligible numerically. The form of the NLL correction suggests that the most natural scale 
is that of the transverse momentum squared of the emitted gluon, see for example [29,30]. 
This confirms the choice of fl2.9j ) as a more natural choice of scale. 

2.3 Numerical results for the linear case 

2.3.1 Different versions of the non-Sudakov form factor 

We will present here numerical results for the solution to the linear CCFM evolution. We 
start with the comparison of the two versions of the non-Sudakov form factor ( |2.6| ) and 
( |2.7|) in the case when the coupling is fixed, and set to be equal to a s = 0.2. In figure [| we 
present the value of the unintegrated gluon distribution A(x,k,p) obtained from CCFM 
as a function of the transverse momentum k, and for the fixed value of the momentum 
p = 10 GeV (left plot) and p = 200 GeV (right plot). Different curves in increasing order 
correspond to increasing values of the rapidity. The initial condition for the equation was 
set to be 

A (0) (k,p) = exp(-k 2 /fi 2 )6(p-k) , (2.10) 

with fi = 1 GeV. It is clear that the 'non-regular' form-factor given by (2.6) gives a 



significant contribution in the low momentum regime due to the fact that k < zq so 
that the expression in the exponent changes sign. This gives a substantial enhancement 
at low values of k. On the other hand it is also clear that both prescriptions for the form- 
factors coincide in the large k regime. One observes that the slope of the k distribution 
becomes steeper when k becomes larger than p, meaning a suppression of the momenta due 
to angular ordering. Momenta larger than p are nevertheless allowed. This change in the 
slope in k is crucial for the analysis with the boundary and the behavior of the saturation 
scale as we shall see later. 

2.3.2 Fixed vs running coupling 

In figure ^ we perform the comparison between the running and fixed coupling scenarios. 
The running coupling is regularized at k 2 . = 0.7 GeV with frozen scenario. As expected 
the running coupling gives the higher contribution in the region of the infrared momenta. 
It falls then below the fixed coupling at high momenta but such that k < p, as the running 
coupling value becomes smaller than the fixed value a s = 0.2. The effect is most pronounced 
for the calculation with p = 200 GeV. We observe though that in the region where k > p, 
the scenario with the running coupling gives results which tend to be larger. One needs to 
take into account the fact that in this region we have stronger relative suppression from the 
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Figure 5: Solution to the CCFM equation (2.2) for fixed a s = 0.2. Curves increasing in magnitude 
are for rapidities Y = 4, 6, 8, 10 respectively. Two values of the external scale are used, p = 10 GeV 
(left plot) and p — 200 GeV (right plot). Dashed blue lines correspond to the non-Sudakov form- 



factor (|2.7|) while solid red curves correspond to the non-Sudakov form-factor (2.6) 
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Figure 6: Solution to the CCFM equation (2.2). Curves increasing in magnitude are for rapidities 
Y = 4,6, 8, 10, 12, 14, 16 respectively. Values of the external scale are p = 10 GeV (left plot) and 
p = 200 GeV (right plot). Dashed blue lines correspond to the fixed coupling case with a s = 0.2 
while solid red curves correspond to the running coupling case. 



non-Sudakov form- factor (we shall discuss this more below). The value of the form- factor 
is smaller for the larger value of a s , which is for the case of the fixed coupling. This is why 
the running coupling calculation tends to be less suppressed in the asymptotic regime of 
the large transverse momenta. 
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Figure 7: Solution to the CCFM equation (2.2). Curves increasing in magnitude are for rapidities 
Y = 4,6,8, 10, 12, 14, 16 respectively. Value of the external scale p = 1 GeV (left plot) and p = 
200 GeV (right plot). Dashed blue lines correspond to the running coupling with scale q; solid red 
curves correspond to the running coupling with scale k. 



2.3.3 Comparison between different versions of the running coupling 

Next let us turn to the implementation of different scales in the running coupling. In figure 
[7| we show the comparison between the standard results and those obtained using q as 
scale. The differences between two choices are in general not too large, with the choice of 
q giving smaller values for the solution. In the case when p S> k there are more and more 
chains with q 3> k contributing, and so typically a s (k) > a s (q). Therefore the solution 
with q as a choice of the scale is slightly lower. 

On the other hand, when p is small (~ 1 GeV) we see that the differences are larger, 
and especially the Y dependence and the slope of A seem to be modified. In this regime 
the phase space is strongly constrained by the maximum angle allowed for the emissions. 
The phase space restriction applies to the real momenta q, while the virtual momenta k 
are not directly restricted. Then in an event where fc>pwe would expect that generally 
k 3> q as well (except early in the chain). This implies that a s (k) < 6t s {q) and one could 
therefore expect to see a faster growth when q is chosen as the scale. However, the opposite 
behavior is observed. This is because a larger coupling also implies that there is a larger 
suppression from A ns . We saw a hint of this behavior already in figure H. When k > p the 
phase space for the real emissions is more strongly constrained, on the other hand ( |2.6j ) 
or (|2.7| ) are independent of p, so even if the real phase space is strongly constrained, the 
form- factor is not. This implies the possibility that a larger coupling in this regime actually 
slows down the evolution due to the large suppression from A ns . 

To further demonstrate that this is indeed the case we show in figure |8| results obtained 
from two different values of the fixed coupling, a s = 0.15 and a s = 0.3, for p = 5 GeV. 
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Figure 8: Solution to the CCFM equation (2.2). Curves increasing in magnitude are for rapidities 
Y = 8,12,16 respectively. The value of the external scale is p = 5 GeV. Dashed blue lines 
correspond to a s = 0.3 while the solid red curves correspond to a s = 0.15. 



For clarity of demonstration we show only solutions for Y = 8, 12, 16. Of course in this 
case since a s = 0.3 > a s = 0.15 for any k, in the region where k < p the solution with 
the larger coupling indeed grows faster (the phase space restriction is not so important 
here). However, we then clearly see that as k > p, the suppression from A ns , which is 
much larger when a s = 0.3, starts to dominate and therefore the solution with the larger 
coupling is strongly suppressed. Consequently in the region k S> p, the solution with the 
larger coupling actually grows slower. 



3. CCFM in the presence of saturation: analytical results 

Having reviewed some of the basic concepts in CCFM we now move on to the question 
of implementing saturation effects and finding the relevant saturation scale Q s . The im- 
plementation of saturation effects will be done using the same approach as in [14, 15], but 
most importantly because in this work we do keep the full dependence of A on the addi- 
tional scale p, we shall discover many interesting new results not noticed before. The main 
implementation is numerical and all the numerical results are to be presented in section H 
In this section we will analytically study the general behavior of Q s . 

To study the saturation scale it is convenient to take the Mellin transform over the 
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momentum k. Following [31] we shall write the Mellin representation of A as 

A(Y,k,p) = J^^)y-(i-i)p H (p/k) . (3.1) 

where Y = In 1/x, and we define p = ln(A; 2 /Qq) with Qq an arbitrary scale. We will 
only consider the fixed coupling case in the analytic formulas. The function H encodes the 
dependence of A on the parameter p. What we know generally is that H ~ 1 when p/k S> 1, 
since the restriction coming from £ becomes irrelevant in that case. This, however, does 
not imply that angular ordering is not important anymore. More precisely it means that 
the ordering in the final step of the evolution is automatic when p> k, but in the previous 
steps of the evolution it is still important. This will reflect itself in the characteristic 
function uj which in the limit p 3> k is given by 

1-7 \ 

w = «• / ^ I ( ^-^ ) H (?/l* + q\)-0(k- q) H(q/k) . (3.2) 




To find the function H one can differentiate the integral equation (|2.2j ) with respect 
to p as done in [31]. One then finds that 

^^ = ^l^ 2 0(Y-ln q /p)e( q -p)A(ln q /p,k, q ) (^ (^J W*0 , 

(3.3) 

where as before k' = \k + q\. Assuming that Y is very large and the solution is in the 
asymptotic region one can neglect the first theta function, which is what is done in [31]. 
However, one should then notice that with the non-Sudakov in ( |2.6| ) which is allowed to 
exceed unity, one would run into a divergence. This is so because when Y — > oo, even if p 
is finite, the phase space for the emissions contributing to the scattering becomes infinitely 
large, and there is a UV divergence appearing when q — > oo. On the other hand such a 



divergence does not appear if the restricted A ns in (|2.7|) is used instead. The reason for 
this is that the choice in ( p.7|) effectively implies that k plays the role of a UV cut (but only 
for the soft emissions). If the restriction from £ is removed (as when Y —> oo) then the 
form factors S e ik and S ne need to be regularized by some UV cut-off. When these factors 
are multiplied together, because of the different signs, they cancel each other in the region 
above k, and thus the dependence on the UV cut-off vanishes. However, the real emissions 
need still be regularized by the UV cut-off, since otherwise the inclusive summation over 
the soft emissions for the given hard emission q would give a divergence as q — > oo. The 
fact that the regions Ai in figures y and |I] are neglected in the derivation of ( [2.7D means 
that one is cutting off all real soft emissions above k which is why k plays the role of the UV 
regulator 6 for the soft emissions. Thus for simplicity we shall use this form factor in the 



6 At this point it should, however, be mentioned that there is no reason for the soft emissions to be 
bounded by k. This would essentially be true if the soft emissions were k± conserving in the angular 
ordered cascade. The fact is that they are k± conserving, but not with respect to the angular ordering, but 
rather to the energy ordering whereby each emission is indexed according to its energy, and the "propagator" 
momenta are defined accordingly. 
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following, though a proper renormalization procedure for all the emissions would probably 
be called for. 

After these remarks, let us now return to the equation for H above. The saturation of 
the solution can be analyzed in the two different regimes of large and small p. 

The case when p > k: 

In case p > k, since q > p (because of the theta function in ( |3.3|) ), we also have that 



q > k. We will then approximate k' rs q. Equation (|3.3f ) is then easily solved to give 

H{p/k) = 1 - ^^ -(1 - e'^-W) (%X^ 

V ' ' (l- 7 )(^ + 2-2 7 ) V >W) 

= i _ ^W (i_ 7 ) W - p) ( 4) 

(l-7)(w + 2-2 7 ) l > 

where in the second equality we neglected the term exponential in Y, because it is very 
small when Y is large. We have also defined p p = Iii(p 2 /Qq). Consistency of the solution 
requires that 

# (!) = 7— ^T ' where /(7) = (1 " 7)(w + 2 - 2 7 ) . (3.5) 

1 + /(7) 

We see that when (p p — p) S> 1, that is when p is very large, H reduces to unity. Since 
(p p — p) is large we can rewrite H as 

H(p/k) = exp f_^B) e -(i-7)(ft>-pA . (3.6) 

The Mellin representation then reads 

^(y, fc >P ) « / — e« y -^)'- i W l *- (1 - r)C *- rt . (3.7) 

7 2vri 

To find the saturation scale, one can essentially follow a line of constant density [16,17] 
which means that the relevant saddle point for the saturation problem is determined by 
the equations (defining p s = ln(Qs/Qo)) 

co s Y - (1 - 7s )p s - ^^ e-V-T^-'J = , (3.8) 

Is 

^Y + Ps-^^e^-^^(^-f s+Pp -p^ =0. (3.9) 

Here all derivatives are with respect to 7 and the second equation is the usual saddle point 
condition while the first equation states that the integrand around the saturation scale is 
constant and of order 1. On the other hand we recall that in the saturation problem for 
BFKL the corresponding saddle point equations read 

co s Y - (1 - ls ) Ps = 0, (3.10) 

J s Y + Ps = 0. (3.11) 
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What we notice, besides the fact that the equations in CCFM look much more compli- 
cated than in BFKL, is that equations ( |3,8|) and (|3.9| ) imply that the saturation anomalous 
dimension 7 S , and thus also the characteristic function uj s , depend generally on p and Y. 
In BFKL on the other hand, 7 S is a pure number and uj s is independent of Y which we 
see from the corresponding equations. As p p — p s — > oo, we see that the structure of the 
CCFM saddle point equations reduce to the BFKL ones, though at finite a s , j s and uj s are 
still different than in BFKL. When a s —> 0, the CCFM solution reduces completely to the 
BFKL one. 

From equation (|3.8|) we can find p s . Let p s be the solution when p p — p s — > oo, that 
is 

(o) 

/4 0) = ^o)^. (3-12) 

1 7« 

Since the general form of uj which depends on p is complicated, we will in the following 
simply make the approximations to set 7 S = 7s and uj s = u s . Then let us write 7 
Ps = Ps - Ps , and we get 

a s H(l) ^_ (1 _ 7s)(pp _ ps) ~ ps = _ (1) 



/, 



e -^-is){ Pp -Ps) e - Ps = 0j where p s = (l- 7s ) pW . (3.13) 



Since (a s H/f)e ^ ^ 3 ^Pp p^ <c 1 (it is suppressed by both a s and the exponential) we 
must have p s <^ 1, and therefore 

p s ~ ^H{1) c -(i-^)( Pv - Ps ) L _ a s H{l) e _ (1 _ 7s)(pp _ ps) \ ^ ^ ^ 

Js \ Js J 



Thus we find that 



w.l 0) .. a„H(i) _m_„(o)n„ ,,.,(o) 



1 - 7i 0) d-7i 0) )/i 0) ' ( } 

which is valid as long as p p — p s 3> 1. The saturation saddle point is determined by 






-4 0)/ = -^ , (3-16) 



as in BFKL. Of course this is the lowest order approximation (lowest order with respect 
to the corrections in p, but it is valid at finite a s ) but since the more general case is 
rather complicated we shall use this formula. In the later section we will determine the full 
saturation momentum via the numerical solution. 



In order to be able to use the formula (3.15) one needs to have a solution for the char 



acteristic function in CCFM. This is essential for obtaining col and 7s . Unfortunately it 
is only possible to extract the characteristic function numerically, as was done in [31]. To 



7 The minus sign in front of the correction term comes from the fact that a finite p reduces the phase 
space and should therefore reduce p s too. 
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Figure 9: Logarithm of the saturation scale from ( 3.15 ) as a function of rapidity. The curves in 
increasing order are for external scales p = 10, 20, 50, 100 GeV. The scale Qq = 1 GeV. The curves 
are truncated at the limit of the applicability of the approximation used. 



estimate here the dependence of the saturation scale on p and Y from ( |3.15| ) we will use a 
simple model for the characteristic function, 



X (7, w) = 2* (1) - *( 7 ) - *(1 - 7 + u) 



(3.17) 



This model corresponds to the LO BFKL equation with the kinematical constraint included. 
It is well known that it leads to the significantly reduced value of the intercept [32,33] and 
thus it should be fairly close numerically to the CCFM characteristic function. One can 



,W 



extract 7^ and ojg solving numerically equation ( |3.16| ). For the standard choice of the 



„(°) 



m 



'(0) 



coupling constant a s = 0.2 we obtained 7s = 0.44, tug = 0.37 and — tOs^' = 0.66. Using 
these extracted values we can determine saturation scale from the equation ( 3.15J ). We 
demonstrate the results in figure ||. The curves increasing are for different values of the 
momentum scale p = 10, 20, 50, 100 GeV. 

It is evident that for small rapidities the growth of the saturation scale is dominated 
by the term exponential in rapidity and the characteristics is similar to the BK equation 
(or BFKL with the saturation boundary). However for fixed value of scale p the depen- 
dence of p s on rapidity becomes nonlinear. The smaller the value of p is, the lower is the 
critical value of rapidity for which the behavior starts to become nonlinear, and there is 
an additional dependence on p. The formula cannot be trusted in the region where p s 
falls off with rapidity. This region is outside the validity of the original approximation in 
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which equation ( 3.15| ) was derived. One would need to compute higher order corrections 



and possibly resum them in order to predict the asymptotic behavior of Q s with Y. 

Thus the leading term for the saturation momentum when p is very large has the same 
structure as in BFKL, with the "speed", — uj' s , of p s being much slower, around 0.65 for 
CCFM as opposed to around 0.98 for BFKL at LO and a s = 0.2. The corrections coming 
from finite p slow down the saturation scale and from ( 3.15J ) we see that the p dependence in 



Q 2 S enters as a power inside an exponential. Summarizing, we observe that the corrections 
from the phase space restriction in CCFM introduce an increasingly complicated shape of 
Q s depending both on p and Y when going beyond the BFKL regime, i.e. p> k. 

Using the results above one can also study the behavior of A near saturation, that is 
when < p — p s < 1. This is done by expanding the function in the exponential of the 
Merlin representation to 0{(p — p s ) 2 ). In BFKL this leads to 

Abfkl e -(i-7«)(A-p.) e - a^gr- . (3.i 8 ) 

The BFKL result features the well known characteristics of the small- x evolution. One is 
the property of geometric scaling which is the fact that the Y and k± dependence enters 
the gluon distribution in the combination k±/Q s , that is via p — p s . Indeed as long as 
p — p s <C yj2uj'lY the dominant dependence of A is on p — p s . A more careful treatment 
performed in [16] modifies this solution so that the Y dependent pre-factor coming from 
the Gaussian integration is replaced by a factor (p — p s + A) for some number A. In this 
case geometric scaling is even more obvious. The second property is the diffusion in k± 
space. This is governed by the Gaussian and the diffusion radius grows as W2u£Y. 

One can perform the exact same manipulations also in (|3.7|) . However the resulting 
expression is rather long and not very transparent so we do not write it down here. Again 
in case p p — p 3> 1 the structure reduces to that in ( |3.18|) . If the p dependence is kept, 
both of the two characteristics mentioned above are modified. For example, diffusion is re- 
duced due to angular ordering, and the diffusion radius obtains a complicated dependence 
on p p — p. Geometric scaling is also not transparent anymore, because of the additional 
p and Y dependence. Besides, generally we would need to take into account the p and 
Y dependence of 7 S and oj s which further breaks the scaling behavior. One may find a 
modified scaling behavior even if it is not as transparent as in the BFKL case, but given 
the complexity of the solutions we shall refrain from pursuing this question analytically. 

The case when p < k: 

The opposite region where p <C k is more complicated, but also turns out to be much 
more interesting. In this case we can no longer use expression ( |3.2| ) for co, and the p 
dependence is stronger. Note that in the argument above, which follows the one in [17], 
saturation is not really explicit. It is understood that the solution represented by the 
Mellin transform is valid only as long as one is above Q s , but the effects of saturation 
are kept implicit. Thus one could actually get the same result if one simply took linear 
BFKL and followed a line of constant amplitude, though of course in the linear case this is 
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rather ad hoc since there is nothing special about any particular line of constant amplitude, 
nevertheless the answer for Q s would be the same. Therefore the leading asymptotic Y 
behavior of Q s extracted in this way from BFKL, without actually explicitly applying 
saturation, is the same as that extracted from the non-linear BK equation. To obtain 
correctly the non-leading Y behavior, however, one has to go beyond this approach but 
what was shown in [16] is that one can get the correct answer from the linear equation if 
the effects of saturation are introduced via the boundary prescription. 

In CCFM the effects of saturation will be seen to be rather dramatic. From the 
numerical analysis we shall see that the phase space restriction coming from p, which 
has its origin in the coherence of the QCD radiation, when combined with saturation, 
which imposes a restriction in the region below Q s , implies that the saturation scale itself 
saturates at a region Q s S> p. This behavior is missed by the naive analysis above since 
one must explicitly apply saturation to possibly see it. A treatment like in [16] may lead to 
an analytic understanding of this novel behavior, but given the complexity of the equations 
involved this may be very difficult to pursue. We shall therefore not apply the analytic 
treatment to this case, even though finding H itself is easy, as done in [31]. Instead we 
move on to the numerical treatment. 

4. Description of the boundary method 

The boundary method applied to the linear BFKL evolution, was shown in [14, 15] to suc- 
cessfully reproduce the results obtained from the non-linear BK equation. The saturation 
momentum Q s and the shape of the gluon distribution for k > Q s were correctly repro- 
duced by this method. The basic idea is to prevent the gluon occupation number from 
growing too large. First, we define a line of constant occupation number via the condition 

A(Y, Pc ,p)=c, (4.1) 

where the number c is of order 1 or smaller (but not much smaller). As obvious from 
the definition, p c depends on both Y and p. In the linear evolution, A will continue to 
grow indefinitely beyond p c when increasing Y. We then enforce a boundary condition 
preventing such an unrestricted growth beyond the critical line p c . In practice this is done 
by applying the boundary for all momenta p < p c (Y,p) — A. Here A is a given number. 
Strictly speaking in CCFM A should not be a pure number but should depend on Y and p. 
However, given that explicit analytic solutions of CCFM are not available, it is not really 
feasible to analytically find the possible dependence of A on Y and p. We will therefore 
for simplicity let A be a constant as in BFKL. The parameters c and A can be viewed as 
free parameters although in BFKL it is known that A ~ ln(l/c). In particular the results 
should not depend on the precise values of these parameters, apart from differences in 
normalization. Our default choice will be to set A = 2 and c = 0.4. In [14, 15] the default 
choice of the boundary condition was to let 

A(Y,p,p) = 0, for all p < p c (Y,p) = p c (Y,p) - A . (4.2) 
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In this case one has a totally absorptive boundary, and this was indeed also the choice in 
the analytical analysis in [16]. This condition is convenient for the analytical analysis in 
the BFKL case because then the problem can be formulated as diffusion in the presence of 
an absorbing wall [16]. 

The precise choice of the boundary condition, should (at least asymptotically) not 
make any difference beyond a normalization factor. Indeed , it was explicitly demonstrated 
in [14, 15] using the numerical analysis that the method is robust for any Y (thus not only 
asymptotic ones), so that different choices still give similar results. In this paper we shall 
not be using the above condition, but we will rather investigate different alternatives. 

One possible choice of boundary is to let A grow logarithmically beyond the line p c . 
Then we let 

A(Y,p,p) = A(Y,p c ,p) + p c (Y,p) - p, for all p < p c (Y,p) . (4.3) 

This is similar to the behavior observed in the BK equation. 

Another choice is to simply freeze the gluon distribution at the critical point, 

A(Y,p,p) = A(Y,p c ,p), for all p < p c (Y,p) . (4.4) 

Once the boundary has been implemented, one can extract the saturation momentum 
from the evolution. Q s can be defined as a line of constant density. One can define p s 
as any line of constant density between p c and p c . The precise choice affects only the 
normalization of Q s . 

We will investigate both choices of the boundary condition. This is because one does 
not really know the details of the non-linear dynamics in the CCFM, and we want to check 
the robustness of the obtained results with respect to the change of the boundary condition. 

5. Numerical results in the presence of the saturation boundary 
5.1 Fixed coupling results 



In figure 1C we show the solution to the CCFM with the saturation boundary implemented 



in the form of freezing as in (kL4|). Here, the coupling is fixed a s = 0.2. We compare 



solutions obtained using (|2.6|) and (2.7) respectively. We see that once the saturation 
effects are included, the differences between the solutions using the different form-factors 
are not that large, as the sensitive infrared region is regularized by the saturation boundary. 
The dependence on p can be examined by comparing the left and right plot in figure [l(]. We 
observe that, similarly to the linear case, the distributions change slopes when k become 
comparable to the external scale p. For k > p, the gluon distribution is effectively cut off 
by the scale p. 

In order to analyze the asymptotic behavior of the solutions we extend the rapidity 



range up to Y = 50. This is shown in figure |ll|, where we use the non-Sudakov form-factor 



from equation ( |2.6[) . We observe the interesting effect that the front of the solution no 
longer travels to higher momenta with uniform 'speed' governed by the saturation scale, 
but that it rather stalls at high rapidities. The origin of this behavior is rooted in the 
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Figure 10: Solution to the CCFM equation (2.2) with saturation boundary. Curves increasing 
in magnitude are for rapidities Y — 4, 6, 8, 10, 12, 14, 16 respectively. Value of the external scale 
p = 10 GeV (left plot) and p = 200 GeV (right plot). Dashed blue lines correspond to the form- 



factor fl2.7|); solid red curves correspond to the form factor (2.6) 
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Figure 11: Solution to the CCFM equation ( J2.2| ) with two different saturation bound- 



aries; (13) (left) and (4.4) (right). Curves increasing in magnitude are for rapidities Y 
10, 15, 20, 25, 30, 35, 40, 45, 50 respectively. The value of the external scale is p = 10 GeV. 



presence of the additional scale p. As the rapidity grows, the solution diffuses to the higher 
values of the transverse momenta k. At the same time the saturated region grows also with 
energy. The unsaturated region is eventually pushed to the region where k 3> p, which is 
then cut off by the angular ordering condition. Therefore the phase space for new emissions 
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Figure 12: The gluon distribution in the presence of the saturation boundary (4.4) for Y 



6,8, 10, 12, 14. The solid lines are obtained using the form factor (2.6) while the dashed lines are 



obtained using (2.7). Left: p = 1 GeV. Right: p = 200 GeV. 



is squeezed, due to the presence of two scales Q s and p. We have done the analysis with two 
different boundaries ([O) (logarithmic growth at low momenta) and ( [4.4[ ) (freezing) and it 
is clear that the front stops to diffuse to higher momenta in both cases. The stalling of the 
front therefore is a robust feature of CCFM with saturation as it does not depend on the 
way the boundary is implemented. We note that the complete stalling of the front (in the 



case presented in figure 11) occurs for momenta k ~ 10 2 GeV, which is order of magnitude 
larger than the p ~ 10. In consequence, the saturation scale Q s starts to behave differently 
with rapidity in this regime. In particular the value of In Q s will no longer scale linearly 
with rapidity but slower than that, and eventually the saturation scale will 'saturate' itself 
at large Y. The complete analysis of the behavior of the saturation scale is presented in 
the later section. 

5.2 Running coupling results 

Let us now turn the case of a running coupling. In figure O we show the different results 



obtained using the two different form factors, (2.6) and ([2,7|). Here, the boundary condition 
Q4.4D is used, and k is used as scale in a s . As in the fixed coupling case we see that the 
differences are much smaller than in the linear evolution, and almost identical results are 
obtained. In particular we have checked that the small difference between the results stays 
the same up to very large values of Y (far beyond phenomenologically relevant values). We 
find that, apart from the region of smallest Y in the figures, the difference is a normalization 
factor independent of p and Y. Consequently the different form factors lead to almost 
identical saturation scales which we shall look more closely at in the next section. 

Since in the running coupling case the evolution is generally slowed down, one needs to 
go to higher values of Y to see the stalling of A. For larger p in fact, the Y values required 
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Figure 13: Solution of A in the asymptotic region of Y between 46 and 70 units. The non-Sudakov 
in (|2.7[) is used and p = 1 GeV. 



for the stalling are so extremely large that we have not reached them in our numerical 
solution. Starting from smaller p values, however, we clearly see how the evolution is 
eventually slowed down to a point where A starts to stall in a region between p c and p c 
(i.e. around p s ), and one can also see how this region then grows with Y. Increasing p we 
see that one needs ever higher Y to see this behavior (as in the fixed coupling case). 

In figure [l^ we show the extension of the solution shown in the left plot of figure 12 
(that is when p = 1 GeV) to asymptotic values of Y between 46 and 70 units. In this 
case only the solution using ( |2.7| ) is presented, the solution using (|2.6| ) looks essentially the 
same. Compared to the fixed coupling case, the region where A stalls is much smaller, and 
it also grows slower. Therefore the saturation of the saturation scale will set in much later 
(to be shown in the next section). To be able to exactly determine at what Y the stalling 
of A starts to become visible, we would need to know the full non-linear formulation of 
CCFM. The application of different boundary conditions seem to, however, indicate that 
this behavior sets in rather late, as manifest in figure [l3|. Moreover there are additional 
effects not taken into account here (such as energy conservation) which will further slow 
down the evolution. 

Let us now turn to the implementation of different scales in a s . In the linear case we 
observed the fact that generally the choice of q gives a slower evolution compared to the 



default choice k. In figure 14 we show the comparison between the results obtained from the 
two different scales. For the larger p values the differences between the two choices of scales 
are very small. For smaller values of p the differences are more pronounced, the solution 
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Figure 14: Comparison of solutions using as a scale in a s : k (solid lines), q (dashed lines) for 
Y = 6, 8, 10, 12, 14. Left: p = 1 GeV. Right: p = 200 GeV. 



with scale q being consistently smaller than with the choice of scale k. As explained earlier 
in section 2.3.3 , this effect is due to the stronger suppression from A ns when q is selected 
as scale. 

This, however, brings us to another point. Namely the fact that one should actually not 
allow any contribution to the non-Sudakov form factor from the region which lies outside 
the region set by the maximal angle £, or equivalently by p. To understand this, remember 
the derivation of A ns shown in figures |3| and f|. The regions Bi in that case were bounded 
by the maximal angle £. When k 3> p, it can be that region Bi partly, or completely, 
lies outside the region bounded by £. In that case also parts of region Cj will be located 
outside the allowed region (however, C% cannot lie completely outside the allowed region). 
In reality, however, one should not include any contribution from those parts of Cj and Bi 
(actually since k > q only figure |3] and the corresponding regions B\ and C\ are relevant). 
Thus one should in equation ( |2.3D actually only include the part of region C\ which is 
bounded by £. By not doing this one is generally overestimating the suppression from A ns 
when k ^> p. The non-Sudakov form factor will then be modified, and it will depend not 
only on k, q and z, but also on p and the previous values of Z{ in the chain. Unfortunately 
this also implies that it is no longer possible to write an integral equation as in (^J) or 
(|2.2| ) because one needs to keep track of all Zi values in each A ns during the evolution (this 
is because A ns then depends not only the vertex k' —> k + q but also on earlier vertices). 
One can only write down the gluon distribution in an unfolded form, explicitly summing 
it over all possible chains. 

Thus when using the scale q instead of k we do find some differences in the evolution. 
On the other hand, for reasonable values of Y these differences are rather small, and for p 
values around 10 GeV and higher they are negligible. We will also see in the next section 
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Figure 15: The saturation momentum Q s for p = 1,10,40 and 200 GeV respectively (from 
bottom to top) as a function of Y, for fixed coupling a s = 0.2 and using the non-Sudakov (2.7). 



Left: Using boundary condition (4.3). Right: Using boundary condition (4.4) 



that the respective saturation momenta are quite similar, despite the noted differences. 
We have also noted that when k > p one should actually modify A ns to be consistent. 
However, the implementation of this modification is probably in practice only feasible in a 
MC simulation where one can keep track of the history of the entire gluon chain. 

5.3 The saturation scale 

Fixed coupling results 



From the solutions to the gluon distribution one can extract the saturation scale Q s . 
We have already commented on the definition of Q s in section [|. The exact value of the 
constant defining Q s only affects its normalization which we cannot control anyhow, the p 
and Y dependences are independent of the value of this constant (of course the constant 
should not be varied over many orders of magnitude). 

In figure 15 we show the results for the two different boundary conditions ( |L3[ ) and 
when the coupling is fixed and the non-Sudakov in fl2.7|) is used. In both cases we 



find a behavior consistent with the expectation in ( 3.1 5| ) which holds as long as Q s < p. 
As Q s grows above p, as is most visible in the case when p = 1 GeV in figure [15|, we see 
that the growth of p s = IiiQ^/Qq on rapidity becomes non-linear. The Y range shown 
in the figure is similar to what is accessed at the LHC. The absolute values of Q s should 
not be considered as realistic estimates of what Q s would be at the LHC (a fixed coupling 
calculation is not so relevant phenomenologically anyway). The two different boundary 
conditions give rise to slightly different Y dependences for Q s for each given p, but the 
differences are rather small. 

The rather dramatic change in the Y dependence of Q s is visible in figure [Eg where we 
extend the results from figure |15| to much larger Y (far beyond experimentally accessible 
values). As we see, the general shape of Q s does not depend on the precise boundary 
being used. We see that once Q s > p, the dependence of p s on Y becomes non-linear, and 
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Figure 16: The saturation scale from figure [fq extended to very large Y in order to clearly see 
the saturation of Q s . The three set of curves are for p = 1, 10 and 200 GeV respectively (from 
bottom to top) and the two different set of points within each set are for the boundary conditions 



(4.3) and (4.4) 



eventually it saturates to some value depending on p. Of course we cannot say exactly 
whether there exist a given constant C such that Q s < C as Y — > oo or whether Q s keeps 
growing, albeit extremely slowly. Nevertheless, it is clear that, the growth of Q s on rapidity 
saturates and that the dependence of p s on Y is highly non-linear. Interestingly, we observe 
that there exists a certain scaling in the sense that once Q s has saturated, the ratio between 
the values for two given values of p is independent of p. In other words, we find that the 
ratio Q s (Y,p)/p is a constant independent of p, once Q s 3> p. The value of this constant 
depends on the precise boundary condition used, for example for the boundary in fl4.3| ) we 
find Q s (Y,p)/p ~ 90 when Y is very large, while for the boundary in ( |4.4| ) we instead get 
Q s (Y,p)/p^40. 

Next we investigate the results obtained using the non-Sudakov in ( |2.6| ). In this case 
since the normalization of A is larger, the normalization of Q s will also be larger. We 
only show the results using the boundary condition ( f4,3| ), the results for ( |4.4j ) are similar. 
The plots for Q s are shown in figure |17] and we find exactly the same type of behavior 
as above. This is further seen from the right figure where the comparison to the above 
results are shown. As can be seen, the energy dependence of Q s is the same regardless of 
the non-Sudakov chosen, and we find the difference in normalization to be a factor roughly 
1.35, regardless of the value of p. 



Running coupling results 

It is well known that the growth of Q s is slowed down by the running of the coupling. 
Physically it is of course consequence of the fact that the evolution leads to the diffusion 
towards higher transverse momenta where the coupling is much smaller. In BFKL, the 
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Figure 17: Left: Q s obtained using the non-Sudakov in (EJJf) and the boundary condition (4.3) 
for p = 1,10,40 and 200 GeV respectively (from bottom to top). Right: Comparison with the 
results from the same boundary condition but using (^Tj) as non-Sudakov for p — 1 and 200 GeV. 



linear dependence p s ~ Ay at fixed coupling is for a running coupling replaced by p s ~ 
y/2X^boY where for N f = 0, y/2\fo « 3.26. In CCFM we again expect to find a BFKL- 
like behavior as long as Q s <}). As a consequence of the much slower growth, one needs 
thus in the running coupling case much higher Y values in order to reach the point where 
Q s > p. This implies that the stalling of the growth, and the consequent saturation of Q s 
is delayed to higher energies. 

We start by again comparing the effects of different boundary conditions. In addition 
to the two boundary conditions presented above, we also try an additional boundary where 
we simply fix A to be a given constant behind the saturation front. In particular we choose 



A(Y,p,p) = l, for all p<p c (Y,p). 



(5.1) 



We will with this choice further demonstrate the robustness of the method. 



In figure |18| we compare different boundary conditions for both of the non-Sudakov 
form factors. We see that the differences are minor and that they can be removed by p and 
Y independent scale factors. For both cases we find a difference in scaling of around a factor 
1.20. We also see that Q s as in BFKL grows much slower and also that, as a consequence, 
the dependence on p is weaker. As in the fixed coupling case, as long as Q s <C p we find a 
BFKL-like behavior, thus in the phenomenologically relevant Y interval shown in the figure 
we find that for p = 200 GeV, p s grows proportional to \Y. As Q s becomes comparable 
to p, however, this dependence is altered and becomes more complicated. Continuing to 
larger Y, one again finds that Q s eventually saturates. From figure ^ we also see that the 
difference between the results for two non-Sudakov form factors are rather minor, and turns 
out to be a pure scaling, around a factor 1.10, independent of p and Y . These results show 
that seemingly very different linear evolutions obtained from the different non-Sudakov 
form factors give almost identical results once saturation is properly implemented. In 
this way any potential ambiguity of the formalism is removed. We note also that it was 
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Figure 18: Q s for different boundary conditions and p values obtained from a running coupling; 



fl4.4p andp = 1 GeV (filled circles), (|y) andp = 1 GeV (open circles), (M) andp = 200 GcV (filled 
squares), (|5.f|) and p = 200 GeV (open squares). Left: non-Sudakov (2.7). Right: non-Sudakov 
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Figure 19: Comparing the results for the two different non-Sudakov form factors. Filled squares 
(p = 200 GeV) and circles (p = 1 GeV) are obtained from (2.6) while the open ones from (2.7). In 



this case the boundary condition is given by (4.4). 



demonstrated in [15] that the explicit inclusion of the kinematical constraint again only 
leads to a different normalization, once saturation is implemented. 



In order to demonstrate the saturation of Q s we extend in figure 20 the running 
coupling results to asymptotic values of Y. Once again we see the complete saturation 
of Q s for smaller p, but in this case one needs much larger values Y to reach the stalling 
region. Furthermore we see that the scaling behavior noted above in the fixed coupling 
case, namely that Q s (Y,p)/p saturates to a number independent than p, is no longer true 
when the coupling is running. Additionally, as in the previous figures, we see that the 
dependence on p is much weaker than in the fixed coupling case. 
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Figure 20: Asymptotic behavior of Q s for a running coupling using the boundary condition (4.4) 



and the non-Sudakov (2.7) for p = 1, 10,40 and 200 GeV (from bottom to top). 



Next, we extract Q s in the case when q is used as scale for the running coupling. As we 
have seen from the results for the gluon distribution, the evolution for smaller p is generally 
slowed down in this case, and consequently we find a slower growing saturation scale. The 
results are shown in figure £l]. Since the evolution is generally slower, the p dependence 
for smaller Y is also weaker. For phenomenologically interesting Y values, however, the 
differences are not that large. In all the results presented here for Q s one should keep in 
mind that the estimates are highly optimistic, in the sense that Q s is rather large already 
at low Y. From phenomenological analysis performed at HERA, we do not expect to find 
such a large normalization of Q s at the LHC. We discuss the phenomenological relevance 
of our method in the next section. 

Finally let us mention that, though not shown, we have also checked that the conclu- 
sions are not altered when a different set of the boundary parameters are used. In addition 
to the default choice A = 2, c = 0.4 we have also used A = 1, c = 1 and A = 5, c = 0.1. 
The different choices lead to different normalizations of Q s but do not alter the Y or p 
dependence. The method is therefore rather general and robust. 

6. Summary and outlook 

In this paper we have analyzed the nonlinear effects in the CCFM equation using a bound- 
ary method. We have found that, the resulting gluon distribution exhibits similar properties 
as in previous BFKL/BK studies only in the region of transverse momenta lower than the 
scale p in CCFM. As the rapidity grows and the momenta diffuse to the higher values, 
the dependence on the rapidity and p becomes highly complicated, and completely novel 
effects are observed. As a result the diffusion of the transverse momenta to higher values 
is significantly suppressed due to angular ordering when the typical scales exceed p. For 
ultrahigh rapidities this results in the stalling of the propagation front of the solution. 
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Figure 21: Left: The saturation momentum when g is used as the scale in a s . The boundary 
condition (4.4) is used, and p = 1,10,40 and 200 GeV (from bottom to top). Right: Comparison 



with the standard choice of k for the scale of a s , for p = 1 and 200 GeV. 



From the gluon distribution we have extracted the saturation scale in CCFM. We have 
found that due to the presence of the additional scale in CCFM, the saturation scale be- 
comes also dependent on this scale. Also, the high rapidity behavior of the saturation scale 
is significantly modified in CCFM. For ultrahigh rapidities the saturation scale stops to 
grow completely, and this is related to the above-mentioned stalling of the gluon distribu- 
tion. The smaller the value of p, the lower are the rapidities for which deviations from the 
BK-type behavior of the saturation scale are observed. Analytical calculations confirm the 
numerical analysis on the behavior of the saturation scale. 

We have investigated the cases both in fixed and running coupling scenarios, and found 
that in the running coupling case the novel features are present, yet for the higher values 
of rapidities. We have also investigated the dependence of the gluon distribution on the 
choice of the non-Sudakov form-factors. In the case when the boundary is present there 
is very little sensitivity of the solutions on the different form factors, due to the fact that 
the nonlinear effects basically regularize the infrared region where the differences are the 
largest. On the other hand, in the linear case we have found that the different forms of the 
form factors lead to large discrepancies of the solutions. 

Let us now comment on the possible phenomenological relevance and importance of the 
method and results presented here. The results obtained here indicate that the saturation 
of Q s occurs at extremely large values of rapidity. In the phenomenologically relevant 
case of the running coupling, the rapidities needed to observe this phenomenon are about 
Y = 30 — 40. Such extremely large values are not accessible at present particle colliders. 
One should also take into account the fact that here we only included the 1/z term of 
the splitting function. The other, nonsingular at z = terms are important too and they 
are known to reduce the speed of the rapidity evolution. One also needs to properly take 
into account the effects of the energy momentum conservation which will further reduce the 
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growth. Taking these effects into consideration we can conclude that it will be very difficult 
to observe the novel effect of saturation of Q s in the near future. One the other hand one 
needs to note that the exact values of the rapidities at which stalling occurs could depend 
on the details of the initial condition. What could be however relevant phenomenologically 
is the dependence of the saturation scale on the additional scale p in the regime where 
Q s is of order p. For example in DIS, in the usual saturation-based approaches one treats 
Q and Q s (x) as completely independent scales. The results of our analysis indicate that 
saturation scale should be correlated with Q when both scales are of the same order. 

The interesting question is, however, whether or not it is urgent to introduce saturation 
effects for studies at LHC involving unintegrated parton densities. This question is relevant 
and important, as saturation affects the shape of the unintegrated distributions, even at 
higher k. Usually the distributions are extracted using linear evolution equations and 
they will be used in almost every physics analysis at the LHC. A recent study [34] has 
reported deviations from small- x resummed DGLAP evolution in the HERA data that 
might be related to saturation effects. Therefore we think that there is a need to analyze 
the possible saturation effects in the evolution and their impact onto the LHC observables 
with substantial precision. 

To study in more detail the relevance of saturation, one would first have to include 
additional corrections to the hard emissions. Obviously the soft emissions associated with 
the 1/(1 — z) pole are important to include as well as the other non-singular terms in z of 
the splitting function. We intend to return to this issue in a future publication. The soft 
emissions are included in the CASCADE event generator, and as already outlined in [14,15], 
the boundary method can be implemented in an event generator much in the same way. A 
first attempt of applying a boundary method in CASCADE was reported in [35]. In this 
case, however, Q s was assumed a priori to be of the standard form Q s ~ x~ x [19], and it 
was then used as a cut-off in k. As we have seen here, Q s in CCFM depends on rapidity 
and p in much more complicated way. Moreover, the power of the boundary method is 
that it gives a Q s consistent with the evolution. 

A problem we see with the application of the boundary method in CASCADE is that 
the results coming from it show almost no growth at all of A from the small- x evolution. 
The relatively steep growth in x needed to fit data at HERA is rather provided by the 
x dependence of the off-shell hard matrix element and the initial condition. Obviously, if 
there is very little growth coming from the perturbative gluon cascade, then saturation will 
not be relevant. 

A similar effect is observed in the resummed BFKL/DGLAP approaches [29,30,36-42]. 
It was consistently found that the resummed gluon-gluon splitting function does indeed 
have the singularity at small x, but the subleading corrections cause the resummed splitting 
function to be very close numerically to NLO DGLAP in the region x > 10 -4 . Therefore in 
the kinematic regime covered by HERA the evolution of the gluon distribution is controlled 
by the DGLAP equations. On the other hand it was found in earlier analysis [32] that 
when evaluating F2 structure function from the kx factorization approach together with 
unintegrated gluon distribution, large enhancement in the low x region is provided by the 
off-shell matrix element, similarly to what was observed in CASCADE. 
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A related problem was also reported in [43] where some phenomenological applications 
of CCFM were presented. In that case, again only the hard emissions were included, but 



to simulate possible non-leading effects, the non-Sudakov (2.7) was modified by simply 
changing the integration intervals, resulting in a larger suppression coming from the form 
factor (the exact motivation for this change was given by checking that the resulting high 
energy saddle point more or less agree with the one extracted from NLL BFKL studies). 
The result was that the small- x growth turned out be much delayed again, and as discussed 
in that paper this delay might be behind the fact that CCFM predictions on forward jets 
seem to underestimate the growth of the jet cross section with decreasing x (for recent 
studies of forward jets within CASCADE see [44,45]). 

On the other hand, a reason why one has to introduce a strong suppression from 
additional effects, which then delay the small- a; growth, is because it is otherwise not 
possible to fit F2 from the linear evolution. This is not so surprising since the linear 
evolution at small momenta gives a very strong growth which we know must be tamed 
by non-linear effects. Even if one only performs a fit to Fi for moderate values of Q, 
because one is integrating A to smaller k, saturation effects are not completely negligible. 
One can therefore imagine that a procedure like that in [43] , where the non- leading effects 
were introduced by hand, easily overestimates the suppression coming from them when 
saturation effects are not properly included. On the other hand the non-leading pieces are 
more properly included in CASCADE which has the same problems as in [43]. An event 
generator is, however, rather complex, including many different components, and it is not 
always very transparent how big the role of the first principle physics in it is. Therefore 
we believe it will be very valuable to include the soft emissions in a simple numerical 
simulation as done here, so that one can more closely study the combined effects of the 
different physics components in a simpler environment. This may shed light on the possible 
relevance of saturation for hard physics at the LHC (obviously for the possible description 
of the Underlying Event (UE) saturation effects such as multiple scatterings etc are very 
important. For a possible description of the very complex UE at the LHC one will need an 
entirely different approach to the problem, for some recent articles on this issue see [46,47]). 
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